In this vignette we will show:
Univariate and multivariate point pattern data analysis methods for imaging-based approaches.
This includes global metrics on the entire field of view and local variants thereof.
The use case is a MERFISH data set from Moffitt et al. (2018).
For more background on point pattern analysis we refer the reader to our more detailed vignettes and Baddeley, Rubak, and Turner (2015).
We assume that cells (or transcripts) can be approximated as points locations.
spatstat package and ppp objectA powerful package to analyse point patterns in R is called spatstat (Baddeley and Turner 2005). The main data object to compute on is called a ppp object. ppp objects describe point patterns in two dimensional space, ppx objects create multidimensional point patterns. A ppp object is made up of three specifications (Baddeley and Turner 2005):
On this central object, various spatstat metrics can be calculated.
SpatialExperiment Object
SpatialExperiment object as introduced by Righelli et al.Often, the starting point in spatial omics data analysis is a SpatialExperiment (or similar) object in R with Bioconductor. The data we consider here is a MERFISH assay of a mouse preoptic hypothalamus (Chen et al. 2015; Moffitt et al. 2018).
suppressPackageStartupMessages({
library(SpatialExperiment)
library(spatstat.geom)
library(spatstat.explore)
library(dplyr)
library(ggplot2)
library(patchwork)
library(reshape2)
library(stringr)
library(tidyr)
library(magrittr)
})# load the data from ExperimentHub
source("../code/load_data.R")
# source some helper functions
source("../code/utils.R")
library('spatialFDA')
theme_set(theme_light())
# load the SpatialExperiment object
spe <- readRDS("../data/spe.rds")
speclass: SpatialExperiment
dim: 161 73655
metadata(0):
assays(1): exprs
rownames(161): Ace2 Adora2a ... Ucn3 Vgf
rowData names(0):
colnames(73655): 6749ccb4-2ed1-4029-968f-820a287f43c8
6cac74bd-4ea7-4701-8701-42563cc65eb8 ...
6b666f81-7b73-4100-9e02-b5381b39f0f3
fdcddd97-7701-462a-b48f-979111245bd5
colData names(7): Animal_ID Animal_sex ... Neuron_cluster_ID cluster_id
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):
spatialCoords names(2) : x y
imgData names(0):We see that we have an object of class SpatialExperiment with \(161\) genes (rows) and \(73655\) cells. This object is very similar to a SingleCellExperiment object except it has the added spatialCoords slot. One slot in the colData is called sample_id which defines the so called z-slices. The three dimensional tissue is cut in the z-axis into consecutive two dimensional slices (Righelli et al. 2022).
Next, we want to extract three slices of this SpatialExperiment object and convert the 2D slices into ppp objects.
# define the Z-stacks that you want to compare
zstack_list <- list("-0.09", "0.01", "0.21")
# small helper function to extract the z-slices and
# convert them to `ppp` objects
selectZstacks <- function(zstack, spe) {
spe[, spe$sample_id == zstack] |>
.ppp(marks = "cluster_id")
}
pp_ls <- lapply(zstack_list, selectZstacks, spe) |>
setNames(zstack_list)
pp_ls$`-0.09`
Marked planar point pattern: 6185 points
Multitype, with levels =
Ambiguous Astrocyte Endothelial Ependymal Excitatory Inhibitory Microglia OD
Immature OD Mature Pericytes
window: rectangle = [1307.7716, 3097.8543] x [2114.368, 3903.386] units
$`0.01`
Marked planar point pattern: 6111 points
Multitype, with levels =
Ambiguous Astrocyte Endothelial Ependymal Excitatory Inhibitory Microglia OD
Immature OD Mature Pericytes
window: rectangle = [1222.5635, 3012.4248] x [-3993.535, -2202.755] units
$`0.21`
Marked planar point pattern: 5578 points
Multitype, with levels =
Ambiguous Astrocyte Endothelial Ependymal Excitatory Inhibitory Microglia OD
Immature OD Mature Pericytes
window: rectangle = [1328.6103, 3124.2492] x [2199.257, 3993.499] unitsWe see that we obtain a list of three ppp objects for the three z-slices \(-0.09, 0.01, 0.21\).
We can plot one of these slices, e.g. slice \(-0.09\) with ggplot
# create a dataframe from the point pattern
pp_df <- pp_ls[["-0.09"]] |> as.data.frame()
# plot with ggplot
ggplot(pp_df, aes(x, y, colour = marks)) +
geom_point(size = 1) +
coord_equal()
As stated above, one important aspect of a point pattern is the observation window. It represents the region in which a pattern is observed or e.g. a survey was conducted (Baddeley, Rubak, and Turner 2015, 85). In most microscopy use cases we encounter window sampling. Window sampling describes the case where we don’t observe the entire point pattern in a window but just a sample (Baddeley, Rubak, and Turner 2015, 143–45).
The window of a point pattern does not need to be rectangular; we can receive round biopsies or calculate convex and non-convex hulls around our sample (Baddeley, Rubak, and Turner 2015, 143–45).
Let’s investigate the observation window for the slice \(-0.09\).
# subset point pattern list
pp_sub <- pp_ls[["-0.09"]]
# base R plot of all marks
pp_sub |> plot()
Here, we have a rectangular window around all points.
Let’s investigate what a round window would look like:
pp_sub_round <- pp_sub
# calculate circle with radius 850 µm and a center at the
# centroid of the window would look like
w <- disc(r = 850, centroid.owin(Window(pp_sub)))
Window(pp_sub_round) <- w
pp_sub_round |> plot()
Correctly assigning windows is very important. The window should represent the space where points are expected. Point pattern methods assume that the process that we observe is a sample from a larger process (this concept is called window sampling). If the window is too small it would lead to a false underestimation of the area where the points can be potentially observed. This problem of where we can observe points and where not (beyond the boundary of the window) leads to a range of problems collectively called edge effects (Baddeley, Rubak, and Turner 2015, 143–45). We will discuss these later.
The next concept that defines a point pattern is that marks can be associated with the points. The points can also have no mark, which we would call an unmarked point pattern
unmark(pp_sub) |> plot()
Marks can be univariate or multivariate variables that are associated with the points (Baddeley, Rubak, and Turner 2015, 147). In the context of cell biology we can distinguish between discrete marks (e.g. cell types) or continuous marks (e.g. gene expression).
In our example, we have a multitype point pattern, meaning there are different cell types that serve as marks for the point pattern. Multitype means that we consider all marks together. The opposite is multivariate, where we consider the marks independently (Baddeley, Rubak, and Turner 2015, 564 ff.).
First the multitype case:
pp_sub |> plot()
Then splitting the point pattern and plotting a multivariate view on the same pattern.
pp_sub |>
split() |>
plot()
Marks can as well be continuous as in the case of gene expression. We choose some genes from the original paper and look at their distribution (Baddeley, Rubak, and Turner 2015, 637; Moffitt et al. 2018).
# subset the original SpatialExperiment to our example slide -0.09
sub <- spe[, spe$sample_id == "-0.09"]
# Genes from Fig. 6 of Moffitt et al. (2018)
genes <- c("Slc18a2", "Esr1", "Pgr")
gex <- assay(sub)[genes, ] |>
t() |>
as.matrix() |>
data.frame() |>
set_rownames(NULL)
# gene expression to marks
marks(pp_sub) <- gexNow that we have points with multivariate continuous marks
# create a dataframe in long format for plotting
pp_df <- pp_sub |>
as.data.frame() |>
pivot_longer(cols = 3:5)
ggplot(pp_df, aes(x, y, colour = log(value + 1))) +
geom_point(size = 0.5) +
facet_wrap(~name) +
coord_equal() +
scale_color_continuous(type = "viridis")
We note that the expression of the genes Pgr and Slc18a2 is very evenly distributed with some elevations in the middle of the structure. Esr1 shows a half-circle like structure in expression. Note that the expression is here log transformed counts offset by one (to avoid problems with log of zero).
We can compare patterns between marks of the same type. This is referred to as a within mark comparison in our vignette. We can compare discrete marks, so the distribution of one single mark, e.g. a cell type.
# create a dataframe from the point pattern
pp_df_discrete <- lapply(zstack_list, function(x) {
df <- pp_ls[[x]] |> as.data.frame()
df$stack <- x
return(df)
}) |> bind_rows()
# select OD Mature cells
pp_df_odmature <- pp_df_discrete[pp_df_discrete$marks == "OD Mature", ]
ggplot(pp_df_odmature, aes(x, y, colour = marks)) +
geom_point(size = 0.5) +
facet_wrap(~stack, scales = "free") +
theme(aspect.ratio = 1)
Here, we plot the distribution of mature oligodendrocytes across three slices of one 3D brain sample.
Continuous marks can be compared as well, e.g. the expression of a gene across slices of a tissue
pp_df <- lapply(zstack_list, function(x) {
# subset the original SpatialExperiment to our example slide -0.09
sub <- spe[, spe$sample_id == x]
# Genes from Fig. 6 of Moffitt et al. (2018)
genes <- c("Slc18a2", "Esr1", "Pgr")
gex <- assay(sub)[genes, ] |>
t() |>
as.matrix() |>
data.frame() |>
set_rownames(NULL)
# gene expression to marks
marks(pp_ls[[x]]) <- gex
df <- pp_ls[[x]] |> as.data.frame()
df$stack <- x
return(df)
}) |> bind_rows()
ggplot(pp_df, aes(x, y, colour = log(Esr1 + 1))) +
geom_point(size = 0.5) +
facet_wrap(~stack, scales = "free") +
theme(aspect.ratio = 1) +
scale_color_continuous(type = "viridis")
We note that the gene Esr1 is differentially distributed across the slices of the 3D sample.
Correlation is a second order quantity that measures the dependence between points (Baddeley, Rubak, and Turner 2015, 199). A famous way to measure this is with Ripley’s \(K\), which is a cumulative function that quantifies the “number of \(r\)-neighbours of a typical random point” (Baddeley, Rubak, and Turner 2015, 204; Ripley 1976).
Global correlation measures quantify the correlation in the entire window. Global Ripley’s \(K\) is defined as:
\[ \hat{K}(r) = \frac{|W|}{n(n-1)}\sum_{i=1}^n\sum_{j=1 \\j \neq i}^n\{d_{ij}\leq r\} e_{ij}(r) \]
In the formula above we note a few things:
The function is normalised by the number of points \(n\) and the window size \(|W|\)
the term \(e_{ij}(r)\) is an edge correction - see the section on border corrections further down in the vignette (Baddeley, Rubak, and Turner 2015, 204).
Ripley’s \(K\) function can be variance stabilised, which is referred to as Besag’s \(L\) (Canete et al. 2022; Besag 1977). The idea behind variance stabilisation is to “uncouple” the relationship between mean and variance. By taking the square root of the function in question, the variance is nearly constant across the function (Bartlett 1947).
\[ L(r) = \sqrt{\frac{K(r)}{\pi}} \]
res_ls <- lapply(list('Kest', 'Lest'), function(fun) {
res <- calcMetricPerFov(
spe,
'OD Mature',
subsetby = 'sample_id',
fun = fun,
marks = 'cluster_id',
rSeq = NULL,
by = c('Animal_ID', 'sample_id')
)
res <- subset(res, sample_id %in% c('-0.09', '0.01', '0.21'))
return(res)
})
p_ls <- lapply(res_ls, function(res) {
plotMetricPerFov(
res,
theo = TRUE,
correction = "iso",
x = "r",
imageId = 'sample_id'
)
})wrap_plots(p_ls, guides = 'collect')
The strongest estimate of association between oligodendrocytes is found for the slices \(0.01\). Slice \(0.21\) does not show such a high degree of association at radii \(\leq300\) as the other two slices. This means that the apparent clustering we see in the distribution of points is mainly due to an overall higher number of cells in slide \(0.21\) and not a higher degree of association per se. The black line indicates the expected \(K\) respectively \(L\) function for a completely spatially random poisson process (Baddeley, Rubak, and Turner 2015, 132 ff.).
Same can be done for the arrangement of Microglia cells across these three sections.
res <- calcMetricPerFov(
spe,
'Microglia',
subsetby = 'sample_id',
fun = 'Lest',
marks = 'cluster_id',
rSeq = NULL,
by = c('Animal_ID', 'sample_id')
)
res <- subset(res, sample_id %in% c('-0.09', '0.01', '0.21'))
plotMetricPerFov(
res,
theo = TRUE,
correction = "iso",
x = "r",
imageId = 'sample_id'
)We note that the Microglia cells are distributed close to a Poisson Process across the different slices.
A similar analysis can be performed for continuous marks. As an exercise, try to calculated a mark weighted correlation function markcorr. You get more information on this function by typing ?markcorr. The mark weighted correlation function is defined as:
\[ k_f(r) = \frac{\mathbb{E}[f(m(u),m(v))|u,v \in X]}{\mathbb{E}[f(M,M')]} \]
where the numerator is the conditional expectation of the marks at location \(u,v\) separated by a radius \(r\) and \(f\) can be any function linking the two marks. The denominator is the expectation of two random marks \(M,M'\) (Baddeley, Rubak, and Turner 2015, 603).
genes <- c("Slc18a2", "Esr1", "Pgr")
gex <- assay(spe)[genes, ] |> t() |> as.matrix() |>
data.frame() |> set_rownames(NULL)
# gene expression to marks
colData(spe) <- colData(spe) |> cbind(gex)
res <- calcMetricPerFov(
spe,
selection = 'Esr1',
subsetby = 'sample_id',
fun = 'markcorr',
marks = 'Esr1',
rSeq = NULL,
by = c('Animal_ID', 'sample_id'),
continuous = TRUE
)
res <- subset(res, sample_id %in% c('-0.09', '0.01', '0.21'))
p <- plotMetricPerFov(
res,
correction = "iso",
theo = TRUE,
x = "r",
imageId = 'sample_id'
)
p
We note that at very small radii the expression of the gene Esr1 shows a high association with itself. Later on, the association is less pronounced than in the slices \(-0.09\) and \(0.01\).
Next to observation window metrics, we can calculate point level statistics as well. One such option is the local indicators of spatial association (LISA). This gives one curve per point in the field of view (Baddeley, Rubak, and Turner 2015, 247–48).
pp <- subset(pp_ls[["0.01"]], marks %in% "OD Mature")
L_odmature_lisa <- localL(pp)
df <- as.data.frame(L_odmature_lisa)
dfm <- reshape2::melt(df, "r")
get_sel <- dfm |>
dplyr::filter(r > 200.5630 & r < 201.4388, variable != "theo") |>
dplyr::mutate(sel = value) |>
dplyr::select(variable, sel)
dfm <- dfm |> left_join(get_sel)
p <- ggplot(dfm, aes(
x = r,
y = value,
group = variable,
colour = sel
)) +
geom_line(linewidth = 1) +
scale_color_continuous(type = "viridis") +
geom_vline(xintercept = 200) +
theme(legend.position = "none") +
ggtitle("LISA curves of slice 0.01")p
These curves could be analysed using tools from functional data analysis such as functional PCA (Baddeley, Rubak, and Turner 2015, 247–48; Ramsay and Silverman 2005).
The same analyses as above can be performed between two cell types. The corresponding functions are called cross functions (Baddeley, Rubak, and Turner 2015, 594 ff.). As an exercise, try to implement (similar to the analyses above) a cross comparison between two cell types of interest. With the provided functions this is possible, just give a function you like as input and a vector with two cell types you wish to compare. You can look at functions in the documentation of spatstat.explore. We provide an example below:
# select OD Mature and Microglia cells
pp_df_odmature_microglia <- pp_df_discrete[pp_df_discrete$marks %in% c("OD Mature", "Microglia"), ]
ggplot(pp_df_odmature_microglia, aes(x, y, colour = marks)) +
geom_point(size = 0.5) +
facet_wrap( ~ stack, scales = "free") +
theme(aspect.ratio = 1)
res <- calcMetricPerFov(
spe,
c("OD Mature", "Microglia"),
subsetby = 'sample_id',
fun = 'Lcross',
marks = 'cluster_id',
rSeq = NULL,
by = c('Animal_ID', 'sample_id')
)
res <- subset(res, sample_id %in% c('-0.09', '0.01', '0.21'))
p <- plotMetricPerFov(
res,
theo = TRUE,
correction = "iso",
x = "r",
imageId = 'sample_id'
)p
We note that there is not a very strong co-localisation indicated by the \(L\) curves between mature oligodendrocytes and microglia cells. If we look at their spatial distribution that makes sense since microglia cells are distributed more or less homogeneously in the respective slices.
Edge effects describe the phenomenon that not the entire point process is observed, but rather only the part within the window \(W\). This means the value of various statistics could be biased along the edges (Baddeley, Rubak, and Turner 2015, 213).
There are many corrections for edge effects that are briefly listed here (Baddeley, Rubak, and Turner 2015, 214–19):
Border correction:
Isotropic correction:
Translation correction:
Point patterns are realisations of a point process. In the analysis we make inferences about the point process.
A point process assumes stochasticity. Therefore, HTS-based approaches are not suitable for point pattern analysis.
There are global metrics for the comparison within a celltype or between celltypes.
There are corresponding metrics for single cells and their interactions.
Point pattern analysis allows for the analysis of continuous gene expression marks as well.
sessionInfo()R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Sonoma 14.7.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/Zurich
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] spatialFDA_0.99.1 dixon_0.0-8
[3] splancs_2.01-44 sp_2.1-1
[5] bluster_1.10.0 spdep_1.2-8
[7] spData_2.3.0 tmap_3.3-4
[9] scater_1.28.0 scran_1.28.2
[11] scuttle_1.10.3 SFEData_1.2.0
[13] SpatialFeatureExperiment_1.2.3 Voyager_1.2.7
[15] rgeoda_0.0.10-4 digest_0.6.33
[17] sf_1.0-14 STexampleData_1.8.0
[19] rlang_1.1.1 spatstat_3.0-6
[21] spatstat.linnet_3.1-1 spatstat.model_3.2-6
[23] rpart_4.1.19 MerfishData_1.2.0
[25] EBImage_4.42.0 ExperimentHub_2.8.1
[27] AnnotationHub_3.8.0 BiocFileCache_2.8.0
[29] dbplyr_2.3.4 magrittr_2.0.3
[31] tidyr_1.3.0 stringr_1.5.0
[33] reshape2_1.4.4 patchwork_1.3.0
[35] ggplot2_3.4.4 dplyr_1.1.3
[37] spatstat.explore_3.2-3 nlme_3.1-162
[39] spatstat.random_3.1-6 spatstat.geom_3.2-5
[41] spatstat.data_3.0-1 SpatialExperiment_1.10.0
[43] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
[45] Biobase_2.60.0 GenomicRanges_1.52.1
[47] GenomeInfoDb_1.36.4 IRanges_2.34.1
[49] S4Vectors_0.38.2 BiocGenerics_0.46.0
[51] MatrixGenerics_1.12.3 matrixStats_1.0.0
loaded via a namespace (and not attached):
[1] spatstat.sparse_3.0-2 bitops_1.0-7
[3] httr_1.4.7 RColorBrewer_1.1-3
[5] tools_4.3.1 utf8_1.2.3
[7] R6_2.5.1 HDF5Array_1.28.1
[9] mgcv_1.9-1 rhdf5filters_1.12.1
[11] withr_2.5.1 gridExtra_2.3
[13] leaflet_2.2.0 leafem_0.2.3
[15] cli_3.6.1 labeling_0.4.3
[17] mvtnorm_1.2-5 proxy_0.4-27
[19] R.utils_2.12.2 dichromat_2.0-0.1
[21] scico_1.5.0 limma_3.56.2
[23] rstudioapi_0.15.0 RSQLite_2.3.1
[25] generics_0.1.3 crosstalk_1.2.0
[27] Matrix_1.5-4.1 ggbeeswarm_0.7.2
[29] fansi_1.0.5 abind_1.4-5
[31] R.methodsS3_1.8.2 terra_1.7-55
[33] lifecycle_1.0.3 yaml_2.3.7
[35] edgeR_3.42.4 rhdf5_2.44.0
[37] tmaptools_3.1-1 grid_4.3.1
[39] blob_1.2.4 promises_1.2.1
[41] dqrng_0.3.1 crayon_1.5.2
[43] lattice_0.21-8 beachmat_2.16.0
[45] KEGGREST_1.40.1 magick_2.8.0
[47] pillar_1.9.0 knitr_1.44
[49] metapod_1.7.0 rjson_0.2.21
[51] boot_1.3-28.1 fda_6.1.8
[53] codetools_0.2-19 wk_0.8.0
[55] glue_1.6.2 vctrs_0.6.4
[57] png_0.1-8 gtable_0.3.4
[59] ks_1.14.2 cachem_1.0.8
[61] xfun_0.40 S4Arrays_1.0.6
[63] mime_0.12 DropletUtils_1.20.0
[65] pracma_2.4.2 fds_1.8
[67] pcaPP_2.0-4 units_0.8-4
[69] statmod_1.5.0 interactiveDisplayBase_1.38.0
[71] ellipsis_0.3.2 bit64_4.0.5
[73] filelock_1.0.2 irlba_2.3.5.1
[75] vipor_0.4.5 KernSmooth_2.23-21
[77] colorspace_2.1-0 DBI_1.1.3
[79] raster_3.6-26 tidyselect_1.2.0
[81] bit_4.0.5 compiler_4.3.1
[83] curl_5.1.0 BiocNeighbors_1.18.0
[85] DelayedArray_0.26.7 scales_1.3.0
[87] classInt_0.4-10 rappdirs_0.3.3
[89] tiff_0.1-12 goftest_1.2-3
[91] rainbow_3.8 fftwtools_0.9-11
[93] spatstat.utils_3.0-3 rmarkdown_2.25
[95] XVector_0.40.0 htmltools_0.5.6.1
[97] pkgconfig_2.0.3 jpeg_0.1-10
[99] base64enc_0.1-3 sparseMatrixStats_1.12.2
[101] fastmap_1.1.1 htmlwidgets_1.6.2
[103] shiny_1.7.5.1 DelayedMatrixStats_1.22.6
[105] farver_2.1.1 jsonlite_1.8.7
[107] mclust_6.1.1 BiocParallel_1.34.2
[109] R.oo_1.25.0 BiocSingular_1.16.0
[111] RCurl_1.98-1.12 GenomeInfoDbData_1.2.10
[113] s2_1.1.4 Rhdf5lib_1.22.1
[115] munsell_0.5.0 Rcpp_1.0.11
[117] ggnewscale_0.4.9 viridis_0.6.4
[119] stringi_1.7.12 leafsync_0.1.0
[121] MASS_7.3-60 zlibbioc_1.46.0
[123] plyr_1.8.9 parallel_4.3.1
[125] ggrepel_0.9.4 deldir_1.0-9
[127] Biostrings_2.68.1 stars_0.6-4
[129] splines_4.3.1 tensor_1.5
[131] locfit_1.5-9.8 igraph_1.5.1
[133] ScaledMatrix_1.8.1 BiocVersion_3.17.1
[135] XML_3.99-0.14 evaluate_0.22
[137] BiocManager_1.30.22 deSolve_1.40
[139] httpuv_1.6.11 purrr_1.0.2
[141] polyclip_1.10-6 rsvd_1.0.5
[143] lwgeom_0.2-13 xtable_1.8-4
[145] e1071_1.7-13 RSpectra_0.16-1
[147] later_1.3.1 viridisLite_0.4.2
[149] class_7.3-22 tibble_3.2.1
[151] memoise_2.0.1 beeswarm_0.4.0
[153] AnnotationDbi_1.62.2 cluster_2.1.4
[155] hdrcde_3.4 BiocStyle_2.28.1 ©2024 The pasta authors. Content is published under Creative Commons CC-BY-4.0 License for the text and GPL-3 License for any code.